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We study the structure and stability of nonlinear impurity modes in the discrete nonlinear 
Schrodinger equation with a single on-site nonlinear impurity emphasizing the effects of interplay 
between discreteness, nonlinearity and disorder. We show how the interaction of a nonlinear local- 
ized mode (a discrete soliton or discrete breather) with a repulsive impurity generates a family of 
stationary states near the impurity site, as well as examine both theoretical and numerical criteria 
for the transition between different localized states via a cascade of bifurcations. 



I. INTRODUCTION 



Nonlinear localized modes in discrete systems, also 
called intrinsic localized modes or discrete breathers, are 
self-localized states that exist due to the interplay be- 
tween lattice coupling and nonlinear effects. Spatially 
localized modes of the discrete nonlinear Schrodinger 
(DNLS) equation, known as discrete solitons, have ap- 
peared in many diverse areas of physics, such as bio- 
physics, nonlinear optics, solid state physics [6] and, more 
recently, in the studies of the Bose-Einstein condensates 
in optical lattices [2] and photonic-crystal waveguides 
and circuits [3]. 

In application to nonlinear guided-wave optics, dis- 
crete solitons have been suggested to exist in nonlin- 
ear waveguide arrays [4], and they have been observed 
a decade later in AlGaAs waveguide structures [5] (see 
also a review paper [6]). Several other issues, such as 
the discreteness-induced soliton dynamics and soliton in- 
teraction, have been investigated. Regarding applica- 
tions, it was recently shown that discrete solitons in two- 
dimensional networks of nonlinear waveguide arrays can 
be used to realize various functional operations, such as 
blocking, routing, logic functions, time gating, etc. [7]. 

On the other hand, localized impurities are known to 
play a crucial role in numerous physical systems, not 
only introducing interesting wave scattering phenomena 
[8], but also creating the possibility for the excitation of 
impurity modes, which are spatially localized oscillatory 
states at the impurity sites [9] . The relevant phenomenol- 
ogy has been recognized as important in a variety of phys- 
ical settings ranging from defect modes in superconduc- 
tors [10] to the dynamics of the electron-phonon interac- 
tions [11] and from the propagation of light in dielectric 
super-lattices with embedded defect layers [12] to defect 
modes in photonic crystals [13]. 

It is crucially important to examine the interplay be- 



tween disorder and nonlinearity. For the continuum prob- 
lems, this interplay is known to lead to the existence of 
symmetric and asymmetric impurity modes and their in- 
teresting stability properties (see Refs. [14,15] and refer- 
ences therein) . Here we concentrate on the study of the 
discrete systems such as the DNLS model. Some earlier 
results for the discrete lattices focused on the scatter- 
ing of localized modes by impurities in the DNLS equa- 
tion [16,17], as well as in the Klein-Gordon model [18]. 
The effect of the impurity on localization properties of 
the nonlinear lattice and stability of localized modes was 
studied, to the best of our knowledge, only for lattices 
with anharmonic coupling [19]. 

However, the earlier theoretical predictions [17] and re- 
cent experimental results [20,21] for the interaction of the 
discrete solitons with structural defects in arrays of Al- 
GaAs optical waveguides suggest that discreteness should 
play a crucial role in the properties and stability of non- 
linear localized modes. In this paper, we study, in the 
framework of the DNLS model, how the interaction of 
a discrete localized mode with a repulsive nonlinear im- 
purity leads to different localized states, as well as ex- 
amine, both analytically and numerically, the stability 
criteria for the bifurcations between different localized 
modes near the impurity site. 

The existence of multiple localized states near the im- 
purity site can be understood by means of simple physics. 
Indeed, in the continuum approximation, a two-hump 
localized mode centered at the defect [22] is known to 
be unstable with respect to an exponential growth of 
antisymmetric linear perturbations [14,15]. The above- 
mentioned instability leads to the motion of the local- 
ized mode away from the defect [15], and for small dis- 
tances £ between the mode center and the defect, the 
effective interaction energy can be presented in the form 
Hi n t — Hq — A^ 2 , where A is defined by the defect pa- 
rameters. However, this situation becomes quite different 
in discrete systems because of the so-called "Peierls re- 



lief" , an effective periodic potential due to the lattice dis- 
creteness. Such potential was first calculated for small- 
amplitude breathers in a nonlinear elastic chain described 
by the discrete nonlinear Klein-Gordon equation [23] ap- 
plying an approach suggested for kinks of the Frenkcl- 
Kontorova model [24]. The potential has the following 
structure, H p ~ H\ — B cos(2tt^ / h) , where h is the lattice 
spacing. Thus, a competition between these two poten- 
tials defines a sequence of stable positions of a localized 
mode near the defect. This qualitative picture is valid 
in the limit of small h and, generally speaking, the prob- 
lem should be solved numerically. Below, we consider 
this problem in the framework of the DNLS model with 
a single on-site nonlinear defect. 

The paper is organized as follows. In Sec. II we in- 
troduce our model and discuss the analytical results ob- 
tained for the family of the stationary states localized 
near the defect site. Then, in Sec. Ill we study numer- 
ically symmetric localized modes and their symmetry- 
breaking instability, including the analysis of the linear 
eigenvalue problem. Section IV includes the results of 
the instability-induced dynamics of localized modes and 
the study of the families of asymmetric localized modes 
which appear as stable states due to a balance of the de- 
fect repulsion and the effective trapping potential of the 
lattice. We also find the families of asymmetric modes 
and demonstrate their stability. Finally, Sec. V con- 
cludes the paper. 



II. MODEL AND ANALYTICAL RESULTS 



We consider a DNLS model with a nonlinear impurity 
on one site of the lattice. To examine the interplay be- 
tween discreteness, nonlinearity and disorder, we study 
the case of a focusing nonlinearity in the waveguide ar- 
ray with a repulsive nonlinear defect. The DNLS equa- 
tion describing our model can be written in the following 
dimcnsionless form, 



C 

lU n + —A 2 U n 



+ \u„\ u n 



(1) 



where A 2 u n = (u n+ i + u„_i — 2u n ), u n describes the 
complex envelope of the electric field in the n-th waveg- 
uide, while the dot stands for the spatial derivative along 
the array propagation direction (which we treat here as 
time t). 

In this study, we select the quintic nonlinearity for the 
defect, in order to differentiate its effect from the lattice 
as well as to model non-Kcrr nonlinear response often 
observed in experiment (see, e.g. Ref. [21]). However, we 
expect that the results will be qualitatively valid for any 
type of the power-law nonlinearity of the order n. The 
generalized nonlinearity of the defect can also appear in 
the theory of photonic circuits [25] . 

The coupling coefficient C can be expressed through 



the lattice spacing h as follows C = l/h 2 , and it charac- 
terizes the effective diffraction in the array. We should 
notice that the lattice spacing h may not coincide with 
the relative distance between the waveguide in an array, 
and these values are connected in a more complicated 
way [26,27]. 



Equation (1) stems from the Hamiltonian 



a . 



(2) 



upon setting its derivative with respect to w* equal to 
iu n . 

The linear waves of the form u n = uq exp(z£li — iknh) 
have the (phonon) spectrum CI = — (2/ti 2 ) sin 2 (fcft./2) 
consisting of a band of negative frequencies (SI < 0). 

In the homogeneous (i.e., no impurity or equivalcntly 
a = 0), continuum (i.e., h — > 0) limit, the model (1) has 
an exact localized solution of the form (in the long-wave 
approximation, i.e. hn <C 1) 



u(x, t) 



rj e 



cosh[?7(x — £)] 



(3) 



where A = w 2 /2, A is the frequency of the breather (and 
should be positive to avoid resonances with the linear 
spectrum of plane waves), and r\ essentially defines an 
effective mass of the breather excitation in the long wave 
approximation according to M e ff — 2rj. 

A natural approach to the study of a discrete model 
is to use the expression (3) as an ansatz solution in the 
Hamiltonian (2) in order to obtain information about the 
static solution of Eq. (1) and its stability. The last term 
of Eq. (2) is then easy to evaluate in the framework of 
the ansatz (3). On the other hand, the first two terms 
(in the parenthesis) of Eq. (2) can be evaluated using 
the Poisson summation formula [281 which states that 
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(4) 



where F is the Fourier transform of /, 

F(k) = -j=J J(x)e ik *dx. 



(5) 



Using the above remarks and formulas for the Hamil- 
tonian (2), we can evaluate H = H(£; {77, h}) (up to con- 
stant, £- independent terms) as the following 



where 



H — Hp + Hi n t, 



4-7T 2 m cos(2m7r£//i) 
h 2 ^ sinh( m-K 2 /rih) 



(6) 
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Notice that the first part of the expression pertaining 
to the summation assumes the correct asymptotic (i.e., 
small h) form used e.g., in Eq. (3.5) in [29], as well as in 
references therein. Finally, it is worth mentioning that 
the series in Eq. (6) has terms which are exponentially 
smaller with respect to the leading order terms and hence 
even keeping a small number of terms in the series should 
yield reasonably accurate results. In the leading order 
approximation, the Hamiltonian can be approximated as 
follows 
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The criterion for obtaining a steady state in the re- 
duced variable formalism of H = H(£;{rj,h}) requires 
that 



dH, 



o. 



(8) 



It is worth noting that the fact that the solutions of Eq. 
(8) will be static solutions of Eq. (1) has been rigorously 
proved recently by Kapitula [30], through the use of a 
Lyapunov-Schmidt reduction. Given that the ansatz of 
Eq. (3) is not exact, this result should be interpreted in 
the orbital sense (i.e., there are orbits in the configura- 
tion space which are close to the approximate solution 
of Eq. (8) in conjunction with Eq. (3), in the appropri- 
ate norm). In particular, in this case one can directly 
observe that £ = will always be a solution of Eq. (8). 
Additional solutions are present close to £ = nh (n £ Z), 
and unstable ones close to £ = (n + 1/2) ft,. In the homo- 
geneous limit, the pulse-like steady states are exactly at 
these points, as is well-known, but now the presence of 
the impurity modifies their exact location, albeit slightly 
due to the exponential decay of the impurity potential 
away from its location at £ = 0. 

Of particular interest is the stability of the solution 
at £ = 0. One naturally expects that when the impu- 
rity is sufficiently strong and repulsive, it will reverse the 
local stability picture close to £ = 0, converting the rele- 
vant local energy minimum into a local energy maximum. 
We should note here that as per the positive definite na- 
ture of the Hessian matrix in the focusing case (for a 
single pulse), the local minima of the energy landscape 
H = -ff (£; {n, h}) will correspond to stable solutions. The 
opposite will be true for local maxima. The Hessian ma- 
trix here is given by OP/dA, where P is the power or 
(squared) L 2 norm of the solution P = \\u\\^ [31]. Recall 
that the L p norm is defined as ||u|| p = (J2 n |u„| p ) 1 / p . 
The L°° norm (used also below) is ||u||oo = max„ \u n \. 
P can be thought of as the number of elementary excita- 
tions present in the field. In the optics literature, this is 
often called the "energy" of the solution. 



Hence, the instability criterion in the reduced, one de- 
gree of freedom dynamical system for the relevant steady 
state will be 



8 2 H 



=o <0, 



(9) 



which can be restated in the form of the asymptotic ap- 
proximation, through the leading order term of the series 
in Eq. (6), as follows a > a c , where 
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rjh 



(10) 



In terms of the full scries, the critical value reads 



8^ 4 



Ctr = 
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which can be easily confirmed to have as leading order 
behavior Eq. (10). 

The result (11) describes the stability threshold of the 
discrete mode localized at the impurity site n = 0, so 
that when the impurity is sufficiently strong and repul- 
sive, it will reverse the local stability of the mode close 
to £ = 0, converting the relevant local energy minimum 
into a local energy maximum. However, the results (8), 
(9) are more general, and they suggest the existence of 
many different stable stationary states located at certain 
distances from the impurity according to the local min- 
ima of the effective energy, the minima are produced by 
the effective periodic potential due to the model discrete- 
ness. For all such states, the physical picture based on 
the reduced Hamiltonian is valid only approximately, and 
more detailed numerical studies are required. 



III. LOCALIZED MODES AND THEIR 
STABILITY 



We now turn to numerical results. We numerically con- 
structed solutions of Eq.(l), for the chain with N = 200 
sites, localized at the impurity site no = 100, using a 
Newton method with an initial condition given by Eq. 
(3), assuming A = 0.5 (n = 1). Upon convergence, 
linear stability analysis of the resulting pulses was per- 
formed. This was done by linearizing using the ansatz: 
u n = u ao i+eexp(iAt)[a n exp(-iut)+b n exp(iuj*t)]; u soi is 
the exact solution (around which we are linearizing) and 
the u)'s denote the linearization eigenfrcquencies. Then 
to O(e) we obtain the relevant (linearization) eigenvalue 
problem 



bt 



= J 



a k 
bt 



where J is the linear stability (Jacobian) matrix of the 
form 



(11) 
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( dFj 8Fj \ 
duj du* \ 
dF* OF* I ! 
duj du* ) 

where F l = -^C(u i+1 +u i _ 1 -2u i ) -ufu* + aufu* 2 8 itTl0 . 
When the corresponding eigenfrequencies u are real, then 
the solution is linearly stable. On the contrary, the pres- 
ence of a mode with a nonzero imaginary part in its eigen- 
frequency denotes the presence of an instability. 
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FIG. 1. Nonlinear localized mode in the model (1) for 
h = 1. Solid line with circles depicts the continuum solution 
(3). Stars and pluses denote the discrete solution (h — 1) at 
a — and a = a c « 0.22, respectively. In the discrete case, 
the solid line is a guide for eye only. 

Figure 1 shows a typical mode profile obtained in differ- 
ent cases. These soliton profiles are similar to those dis- 
cussed in the long-wave approximation in Refs. [14,15,22] . 
For a = (homogeneous case), the circles denote the dis- 
cretization of the continuum profile in accordance with 
Eq. (3). The corresponding exact discrete homogeneous 
solution is shown by the stars. These results are shown 
for the lattice spacing parameter h = 1. Notice, that 
discreteness induces a narrowing and increase of ampli- 
tude with respect to the continuum pulse. Finally, the 
pluses indicate the solution at the onset of the instability 
(that occurs as a is increased), for a = 0.22. An interest- 
ing observation is that at the onset of the instability, the 
pulse has essentially regained its continuum shape but for 
a significantly larger amplitude at the central site. 

A more detailed study of the mode instability that oc- 
curs as a is increased is presented in Figs. 2 and 3. Figure 
2 shows three different examples of the discrete localized 
mode and the corresponding linear stability results. The 
top panel is for a = 0.1, in the subcritical case, where the 
stability of the mode can be determined by the absence 
of imaginary eigenfrequencies (see the top right panel in 
Fig. 2). The middle row shows a slightly supercritical 
value of the impurity strength (a = 0.3, while a c ~ 0.22), 
while a strongly unstable, supercritical mode for a = 0.5 



is shown in the bottom panel. It is important to remark 
here that a perhaps counter-intuitive result of our nu- 
meriCf d investigations is that the instability point is not 
coincident with the point (in parameter space) where the 
symmetric mode becomes two-humped. We have gener- 
ically observed that the mode becomes two-humped for 
considerably larger (than the instability threshold) values 
of the defect strength a. 
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FIG. 2. Different localized modes (left panels) and their 
stability properties (right panels) in the model (1) for h = 1. 
Results of the linear stability analysis are given in the form 
of the spectral plane (ui r ,LOi) showing the real and imaginary 
part of the eigenfrequency, respectively. The top row is for 
a = 0, the middle is for a — 0.3 > a c , while the bottom is for 
a = 0.5. 
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FIG. 3. Threshold of the instability onset for the nonlinear 
impurity mode centered at the impurity site, shown as the 
critical impurity strength a c vs. the lattice spacing h. Solid 
line-the full numerical result, dashed line is given by Eq. (10), 
the dash-dotted is given by Eq. (11). 
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Similar calculations were performed for different val- 
ues of the discreteness parameter h and a two-parameter 
stability diagram is shown in Fig. 3. Above the differ- 
ent curves the (centered at £ = 0) solutions are unstable, 
while the opposite is true below the curves. The solid 
line denotes the exact numerical result, the dashed line 
is the prediction of Eq. (10), while the dash-dotted line 
is the prediction of Eq. (11) including m = 50 terms in 
the summation. One can observe that even though the 
asymptotic result gradually fails (as it should) for larger 
h, the full series prediction of Eq. (11) remains very 
close to the fully numerical result. The difference can be 
well accounted for by the approximations involved in the 
continuum ansatz. 



cies). Perhaps it is even more interesting to examine the 
occurrence of this instability in the light of the top panels 
of the figure. The latter show the L 2 and L°° norms of 
the solution as a function of a. It is noteworthy that the 
occurrence of the instability coincides with the change 
of concavity of the L 2 norm, while it also coincides with 
a maximum of the L°° norm. The above results sug- 
gest two additional numerical criteria that can be used 
to identify the appearance of the instability, namely 
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FIG. 4. Top: the I? (solid line) and L°° (dashed line) 
norms of the localized mode vs. a. Bottom: the relevant point 
spectrum eigenvalues of the problem continued over a. The 
solid line shows the pinning (or translational) mode that be- 
comes unstable (dotted line) after a c « 0.22. The dashed line 
shows the breathing (or edge) mode, while the dash-dotted 
line indicates the band edge of the continuous spectrum. 
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(13) 



IV. 



DYNAMICAL EFFECTS AND 
ASYMMETRIC MODES 



Additional insight on the appearance of the instability 
is given in Fig. 4, where we examine the mode stability 
for different values of a (at h = 1). The bottom panel of 
the figure shows the trajectory of the main eigenmodes 
of the point spectrum of the linearization around the so- 
lution. The solid line shows the translational or pinning 
(antisymmetric) mode, while the dashed line shows the 
edge or breathing (spatially symmetric) mode [32]. No- 
tice that for h = 1 (cf. Fig. 1 of Rcf. [32]), the pinning 
mode has a larger frequency than the breathing mode. 
The band edge of the continuous spectrum which lies at 
A = 0.5 [32] is shown by a dash-dotted line. As a is in- 
creased, it is clear that the breathing mode is "repelled" 
by the presence of the impurity towards the phonon band 
edge. On the contrary, the antisymmetric mode gradu- 
ally approaches the origin and becomes unstable (dotted 
part of the relevant curve) for a > 0.22 (thereafter, for 
a > a c there appears an imaginary pair of eigenfrequen- 



To examine the dynamical development of the instabil- 
ity, once the solution becomes unstable for the supercriti- 
cal values of a > a c , we have performed direct numerical 
simulations of Eq. (1) with an initial condition consist- 
ing of an exact, unstable localized mode for a = 0.3, 
perturbed by small noise (of uniform distribution, and 
amplitude ranging from 0.0001 to 0.1). The scenario that 
is described was found to be generic for such perturba- 
tions of the unstable solution. It was thus found that 
as is shown in the top left panel of Fig. 5 for t = 200, 
the unstable mode sustains a symmetry breaking, which 
cleaves it into another stable localized mode centered es- 
sentially at the neighboring (previous) lattice site ( see 
also the bottom right panel), and to a small- amplitude 
mobile pulse that propagates along the chain and car- 
ries an energy excess. It is worth noting that the original 
(unstable) eigenmode had 3 "main" sites oscillating at an 
amplitude of w 0.9, while the resulting breather is very 
strongly localized at a single site (of amplitude w 1.7). 
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FIG. 5. Development of the mode instability for a = 0.3 
and h = 1. The top left panel shows the mode spatial pro- 
file at t — 200. The top right shows the modulus squared 
of the central site and of the site nearest to it, indicating re- 
laxation to an exponentially localized configuration centered 
around the neighboring minimum (at n ~ 99). The bottom 
left panel shows a space time plot of the evolution of the mod- 
ulus squared of the field, while the bottom right indicates the 
relaxation of the localized mode center from site x c ~ 100 to 
the one with x c « 99. 



The trapping of the unstable mode by the neighboring 
site is consistent with the physical picture described by 
the effective Hamiltonian (6), when the mode overcomes 
the neighboring Peierls-Nabarro barrier and becomes lo- 
calized by the next minimum. However, for smaller val- 
ues of the lattice spacing h, when the Peierls-Nabarro 
barrier becomes negligible, the unstable localized mode 
can start its motion through the lattice under the action 
of the initial perturbation and as a result of the develop- 
ment of the symmetry-breaking instability. An example 
of the latter type is shown in Fig. 6 for the defect strength 
a = 1 and h = 0.5. 




t/10 

FIG. 6. Same as in Fig. 5 but for a = 1 and h = 0.5. 
Instability initiates the motion of a localized mode along the 
chain. 
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FIG. 7. Bifurcation diagram of the asymmetric localized 
impurity modes. The thickest solid line (branch a) corre- 
sponds to the mode centered at no = 100, while the corre- 
sponding dashed line (branch b), to the mode centered at 
no = 99.5. The medium thickness solid line and the cor- 
responding dashed one (branches c and d) are centered, for 
a = at no = 99 and n = 98.5 respectively. The correspond- = 
ing saddle-node bifurcation resulting at the disappearance of 
these branches occurs for a ~ 0.00147. Finally, the thin solid 
and dashed lines (branches e and f) correspond to the cases 
of no — 98 and no = 97.5, respectively. The saddle-node 
bifurcation occurs here at a w 0.0105. The dotted line at 
a — 0.0005 indicates the defect strength for which the mode = = 
profiles of the different branches are shown in Fig. 8. In the 
homogeneous case (a = 0) the array of minima and maxima 
of the Peierls barrier occurs for x c = m and x c = m + 1/2 
(rn e Z). Hence, the points with x c = 100,99,98,97... 
have the same norm (the same is true for the branches with 
x c = 99.5, 98.5, 97.5 . . .) for a = 0. In the presence of the im- = = 
purity (e.g., for a =fc 0), this shift invariance is broken. Thus 
we use the "original" (e.g., for a — 0) center position of the 
branch to classify them when 

One of the most interesting findings of our study of the 
time evolution as a result of the mode instability is the 
effective "relaxation" of the unstable mode towards the 
mode centered at the neighboring (to the defect) site. 
The latter result prompted us to explore the modes in 
the vicinity of the defect. Due to symmetry, our study 
was restricted to modes on one side of the impurity site. 
The bifurcation diagram of the relevant branches of so- 
lutions is shown in Fig. 7 obtained for h = 0.5. We 
have observed the relevant phenomenology to be generic, 
but report it for this value of the lattice spacing that ne- 
cessitates narrower ranges of parameter sweeps (i.e., the 
same observations can be obtained for larger h, but due 
to the weaker coupling, they occur for considerably larger 
values of the defect strength a). 

The natural starting point for our analysis is the con- 
tinuation from the case of an homogeneous lattice (a = 
0), where the stationary modes are well-known to exist 
on a lattice site and between two consecutive lattice sites. 
Hence, in the homogeneous case we consider the branches 
centered at no = 100 (branch a), no = 99.5 (branch b), 
n = 99 (c), n = 98.5 (d), n = 97.5 (c) and n = 97 (f). 
We observe that at the instability point a c w 0.00252, the 
stable branch a and the unstable branch b merge. The 
resulting branch is always unstable thereafter. For any 
other branch apart from the central one, we have found 
that the relevant stable (node) and unstable (saddle) so- 
lutions exist for a parameter interval that depends on the 
distance of the central site of the branch from the defect 
site; naturally, the branches whose central site is more 
remote from the defect arc less affected by it (and exist 
for larger intervals of defect strengths). These branches 
eventually disappear in saddle-node bifurcations, as is ob- 
served for branches c and d and e and f in Fig. 7. A set 
of spatial profiles corresponding to the different branches 
are shown in Fig. 8 for the different branches of Fig. 7 



and for a = 0.0005 (corresponding to the vertical dotted 
line in Fig. 7). 
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FIG. 8. Spatial profiles of the localized modes correspond- 
ing to the branches (a) to (f ) in Fig. 8 for the defect strength 
a = 0.0005. A repulsive impurity is located at the site 
n = 100. 

Disappearance of the relevant branches of static so- 
lutions is consonant with the mobility of the unstable 
modes for smaller values of h observed in Fig. 6 above. 
For smaller values of h, the neighboring branches disap- 
pear for smaller defect strengths, hence allowing for no 
static solutions in the vicinity (in configuration space) of 
the unstable solution, thereby resulting in the mobility of 
the latter (i.e., the solution cannot be trapped in neigh- 
boring potential wells as these have disappeared). How- 
ever, in the opposite case of large h, the lattice is strongly 
discrete and an adiabatic growth of the defect strength 
a leads to the mode motion from one stable state to its 
neighboring state (via an intermediate unstable mode) 
through a cascade of the bifurcation points. 



V. CONCLUSIONS 

We have studied the existence and stability of nonlin- 
ear localized modes in the framework of the DNLS model. 
We have used a variational formalism to obtain the an- 
alytical result for the occurrence of an instability due to 
the interplay between attractive nonlinearity and repul- 
sive impurity. We have examined the variational predic- 
tion by means of numerical bifurcation theory, linear sta- 
bility analysis, as well as direct integration of the discrete 
lattice model. We have found the variational prediction 
(and its refinements) to be in reasonable agreement with 
the numerical results, attributing the discrepancies to the 
variation of the original profile with respect to the pre- 
selected ansatz. We have observed that the instability 
develops much before the solution becomes two-humped, 
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and we have developed numerically motivated criteria for 
the identification of the instability. The study of the dy- 
namical evolution of the unstable modes revealed the po- 
tential for mobile localized modes for smaller values of 
the lattice spacing, while it results in relaxation to the 
neighboring potential wells and the corresponding sta- 
ble, static solutions. Finally, the bifurcation diagram has 
been explored as a function of the defect strength, using 
continuation from the homogeneous lattice case. We have 
found that for the branch corresponding to the defect 
site, the stable and nearest unstable solution branches 
merge at the instability point (producing a thereafter al- 
ways unstable branch) , while the neighboring sites corre- 
sponding stable and unstable branches disappear through 
collision via saddle-node bifurcations. Such a bifurcation 
scenario is in a perfect agreement with the concept of the 
Peierls potential that affects the mode motion in discrete 
systems. 
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